A Percolation Model of Diagenesis 
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The restructuring process of diagenesis in the sedimentary rocks is studied using a percolation 
type model. The cementation and dissolution processes are modeled by the culling of occupied sites 
in rarefied and growth of vacant sites in dense environments. Starting from sub-critical states of 
ordinary percolation the system evolves under the diagenetic rules to critical percolation configura- 
tions. Our numerical simulation results in two dimensions indicate that the stable configuration has 
the same critical behaviour as the ordinary percolation. 
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I. INTRODUCTION 

Rocks in general, particularly sedimentary rocks e.g. 
sandstones, limestones etc., have porous structures. Typ- 
ically such a pore space is a highly branched and inter- 
connected network. Study of the pore structure of sedi- 
mentary rocks is important from a practical point of view, 
in problems such as oil-exploration, ground water flow, 
spread of pollutants etc. 

An interesting property of these rocks is that they 
appear not to have a finite percolation threshold [Q. 
These rock materials are conducting when the pore space 
is filled with saline water. It has been observed that 
these rock samples show finite conductivity even when 
the porosity is less than one percent. This implies that 
a connected network of pores exists in the macroscopic 
length scale, even when the porosity i.e., the volume frac- 
tion of the void space is very little. 

Several empirical laws reflect this property. Archie's 
law Q| connects the conductivity <r(<p) and the porosity 
4> in the following way: 

<x{4>)/a w = a<j) z (1) 

Here, a w is the conductivity of water, a ~ 1 is an empiri- 
cal parameter and z ~ 2 is a non-universal exponent that 
depends on characteristics of the rock structure. This law 
suggests that a finite conductivity exists even in the limit 
of <fi — > and therefore the percolation threshold is zero. 

The permeability K(<fi) of the rock structure is related 
to the porosity <f> through a similar power law, known as 
Kozeny equation Q: 

K^)=c^'/S 2 (2) 

where, z' ~ 3, S a is the specific surface area and c is 
an empirical constant. This equation also suggests the 
global connectivity of the pore space is maintained in 
the <fi — ► limit. 

A physical process which is responsible for achieving a 
connected pore structure at very low porosities is known 
as "diagenesis". Diagenesis is a complex restructuring 



process by which granular systems evolve in geological 
time scales from unconsolidated, high-porosity packings 
toward more consolidated, less porous structures. Forma- 
tion of sedimentary rocks starts with deposition of sand 
grains under water or in air Initially this gives 

an unconsolidated and highly porous ~ 40 — 50% sed- 
iment. Sedimentation is followed by compaction under 
pressure and diagenesis, before the consolidated sand- 
stone is formed from the loosely packed sediment . Di- 
agenesis may reduce porosity by an order of magnitude 
and permeability by as much as four orders of magnitude 
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The final characteristic of the pore network depends 
strongly on the diagenetic process. Sandstones are usu- 
ally formed under water, which contains dissolved salts. 
Depending on the nature of the pore-filling fluids, salts 
may be deposited as crystallites in the crevices or along 
walls of the rock structure, a process called "cementa- 
tion" . Otherwise, portions of the existing solid structure 
may get eroded or dissolved out in a "dissolution" pro- 
cess. The former decreases the porosity of the rock while 
the latter increases porosity. The two processes may take 
place simultaneously. The details of the chemical nature 
of the solid and pore filling fluid determines whether di- 
agenesis leads finally to a stable structure, or to a con- 
tinuously developing structure eventually giving rise to 
caverns of macroscopic size. 

Sahimi had classified the theoretical studies of model- 
ing diagenesis in two ways ra]. The approach of "chem- 
ical modeling" relies on solving continuum equations of 
transport and reactions ignoring the morphology of the 
pore space. The second approach is "geometrical model- 
ing" in which the reaction kinetics and mass transfer are 
ignored. These models start with geometrical descrip- 
tions of initial unconsolidated pore space which evolves 
under simple rules leading to reduction of porosities but 
maintaining the connectivity. For example the model of 
Wong et. al. || starts with a regular lattice in which 
each bond is a fluid filled cylindrical tube of uniform 
radius and conductivity. This system evolves to a ran- 
dom resistor network through a random bond-shrinkage 
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mechanism where randomly selected bonds of the net- 
work shrinks its radius by a constant factor. This model 
maintains global connectivity even in the limit of (f> — > 
and reproduces power law behaviour as in Archie's law. 
A second model of Roberts and Schwartz Q starts with 
a Bernal distribution of dense random spheres of equal 
radii modeling grains. These spheres grow in unison and 
the pore space, i.e. the space not covered by the spheres 
shrinks its volume. This model gives a low but non-zero 
percolation threshold cf> c « 3.5%. The bimodal ballistic 
deposition model (BBDM) jlO| tries to represent the de- 
position realistically, but does not address the problem 
of diagenesis. 
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FIG. 1. (a) On an initial configuration as in (i) if the 
central site is first updated one gets the SC in (ii). However 
if the sequential updating rule (I) (see text) is used one gets 
the SC as in (iii). (b) In parallel updating this configuration 
is locked in a period of cycle 2. 

In our model, we do not take into account the effect of 
chemical reactions explicitly, so this is also a geomet- 
rical modeling of diagenesis. We try to simulate the 
restructuring as it may actually occur in porous rocks 
due to fluid flow. Isolated projected grains on a wall are 
smoothed out modeling dissolution, and a gap or cul-de- 
sac in a solid is hlled by deposition modeling cementation. 
The restructuring involves two processes. Growth of the 
solid phase at sites with a relatively larger number of oc- 
cupied nearest neighbours, to represent cementation, and 
removal or culling of occupied sites which are isolated, or 
have too few nearest neighbours, to represent dissolution. 
This algorithm is a stabilizing process leading to a stable 
structure after several time steps. It may be regarded as 
a self-organizing process as discussed recently by several 
authors . 

Before proposing our model, we briefly describe two 
other physical situations which are closely related to our 
study on diagenesis. In the Bootstrap percolation model 
(BPM) occupied sites of a randomly occupied regular 
lattice with certain probability having fewer than certain 
number of occupied neighbours are successively removed 
O Il8| . On repeated application of this process the sys- 



tem reaches a stable configuration where no further sites 
can be culled. The threshold value of the probability 
at which the stable configuration is percolating is calcu- 
lated. 



FIG. 2. A stable configuration of the diagenetic percola- 
tion for a system of size L = 80 and with m=2. Sites on the 
"infinite" incipient cluster are joined by lines and sites on the 
isolated clusters are shown by filled circles. 



Secondly, consider the nearest neighbour Ising model 
at the zero temperature with Glauber spin-flip dynamics 
in the absence of an external magnetic field. Here the 
direction of a spin follows the direction of the majority 
of the neighbouring spins. In the case when there are 
equal number of up and down spins in the neighbour- 
hood, a spin decides its direction with equal probability. 
Recently it has been observed in that starting from 
an arbitrary random initial configuration of spins this 
system does not reach the global ground states where all 
spins are either up or down but arrive at a frozen two- 
stripe state in a finite fraction of cases [Jl9| . 

In next section we describe our model and also the up- 
dating rules used. Section III describes our results and 
we summarize in section IV. 



II. MODEL 

In our model, the sites of a regular lattice are randomly 
occupied (sj = 1) with a probability p representing pores 
and are kept vacant (sj = 0) with a probability 1 — p 
representing solid grains. This configuration therefore 
models the initial unconsolidated porous structure with 
porosity <f>(p) = p. 
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FIG. 3. The porosity 4>{p) as a function of the initial occu- 
pation probability p for a system size L — 64. The continuous 
curve is a fit to the data having the form given in Eqn. 1. 



0.501 
0.500 
0.499 
^ 0.498 
ft 0.497 
0.496 
0.495 



0.494 

0.00 0.01 



0.02 



0.03 

- -1/v 



0.04 0.05 



o - 
0.06 



FIG. 4. The plot of the percolation thresholds p2 C (L) for 
m — 2 for finite system sizes L as a function of L^ 1 ^ . Using 
v = 4/3, the correlation length exponent for ordinary percola- 
tion, we get the linear fit for large L values. The extrapolated 
value for p2 C is 0.5005. 



The occupation status of a site i depends on its neigh- 
bour number i.e., the number of occupied neighbours 
iii = SjS* where, is the occupation of the j-th neigh- 
bour of the site i. All sites of the lattice are sequentially 
updated according to the following diagenetic conditions: 
(i) Culling condition: Occupied sites having fewer than 
to occupied nearest neighbours are vacated i.e., Sj — ► if 
ni < to (ii) the sites with exactly to occupied neighbours 
remain unaltered i.e., Sj — > Sj if ni = to and (iii) Grow- 
ing condition: Vacant sites having more than to occupied 
nearest neighbours are occupied i.e., — * 1 if > to. 

Starting from the initial configuration the system 
evolves in different time steps following these rules. One 
time step consists of update attempts of all the lattice 
sites. One sweep of the lattice results in another occu- 
pied configuration which is again updated by the same 
rules. This process is continued till the system reaches a 
stable configuration (SC) where no further site changes 
its occupied or vacant status. In general the SC may have 
many clusters of occupied sites. However, there exists a 
percolation threshold p mc of p depending on the value 
of to so that the SC must have a spanning ("infinite") 
cluster of occupied sites for p > p mc in an infinitely large 
system. 

Like the cellular automata models, the updating of the 
sites is important in our problem. Three possible sequen- 
tial updating procedures are as follows: (I) Sites are la- 
beled from 1 to I? from left to right along a row and 
from the first row to the last row. (II) Only sites with 
n- L 7^ to are randomly selected and updated. (Ill) The 
lattice is divided into odd and even sub-lattices and are 
updated alternately but sequentially as in (I). 



For BPM, it has been shown in Jf| that the SC is 
independent of the updating sequence. In contrast here 
the SC does depend on the updating sequence because 
culling at one site may inhibit growth at a neigbouring 
site and vice versa. This is seen by considering the neigh- 
bour numbers at all sites of the lattice. When rn < to 
the culling of the site i reduces the neighbour numbers 
at all neighbouring sites by one i.e., rij rij — 1 where 
as the growth at i enhances the neighbour numbers at 
all neighbouring sites i.e., rij — > rij + 1. Therefore the 
culling at one site may suppress the growth at a neigh- 
bouring site and vice versa. In Fig. 1(a) we show an 
example where two different updating sequences lead to 
different SCs. On the other hand, in a fully parallel up- 
date all sites of the lattice are updated simultaneously 
at a certain time depending on the configuration at the 
previous time. There may arise some situations as shown 
in Fig. 1(b) where a particular cluster of sites never goes 
to a stable configuration but takes two different configu- 
rations alternately in a two cycle periodic state. 

A cluster of occupied sites, in which every site has at 
least to occupied neighbours, is called an m-cluster [ jl4) . 
Imposition of our diagenetic rules imply that the sur- 
viving clusters in SC must be m clusters. An isolated 
cluster of occupied sites in a d-dimcnsional hypercubic 
lattice always has some convex corners on the surface 
with d neighbours. Therefore in the case when to > d+ 1 
these sites are always unstable and therefore, the SC can- 
not have any finite cluster and has only one infinite m- 
cluster. A SC at the percolation threshold for to = 2 is 
shown in Fig. 2, the smallest isolated clusters being of 
size 4. 
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FIG. 5. The average mass of the infinite cluster (cir- 
cle) and the largest cluster 25^ (square) at the diagenetic 
percolation threshold pi c (V) are plotted with the system size 
L. Average fractal dimension is obtained from the linear fits 
shown by continuous lines. The inset shows the variations of 
the local slopes d/(L) and we conclude a value of the fractal 
dimension d f = 1.89 ± 0.02. 



III. RESULTS 

The average fraction 4>(p) of the occupied sites in the 
SC is denned as the porosity of this model. We measure 
this porosity as a function of the probability p and this 
variation is plotted in Fig. 3. We find no trace of any 
system size dependence on this variation. A functional 
form like 

0(p) = l/[l + exp((l/2-p)/Ap)] (3) 

fits very well to this data with a value of Ap = 0.072. 
The data as well as the fit are very well consistent to 
0(1/2) = 1/2 as expected from the symmetry of occupied 
and vacant sites. Compared to the porosity <j){p) = p in 
the initial random distribution of occupied and vacant 
sites, the porosity in SC is reduced by an order of magni- 
tude when p < 0.35. We consider this as the reflection of 
the diagenesis process in nature observed in our model. 

Like BPM, the culling condition in our model does 
not contribute to change the percolation threshold for 
m < 2. For example, nothing is culled for m—0, iso- 
lated sites are culled for m — 1 and the dangling chain 
of sites are culled for m = 2. Since the connectivity of 
the system is not affected by these culling processes, the 
difference between p mc and p c {ord) is due to the grow- 
ing condition for m < 2, where p c {ord) is the ordinary 
percolation threshold. For p values very close to p c (ord) 



but smaller than it, there may be some initial configura- 
tions which are not connected because of the presence of 
only few vacant sites. If these sites have more than m 
occupied neighbours they will now be occupied ensuring 
the global connectivity of the system. Therefore it is ex- 
pected that the p mc < p c (ord) for m < 2 on an arbitrary 
lattice. Therefore as the growth rule helps in attaining 
a connectivity in the system we expect p mc < p c (BPM) 
for any arbitrary lattice and for any arbitrary value of m. 
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FIG. 6. The distribution of the cluster sizes of the stable 
configuration at the percolation threshold of a lattice of size 
L = 2048 and with m = 2. 



The percolation threshold p mc is the minimum value 
of the probability p beyond which an infinite cluster of 
occupied sites exists with probability one in the SC on 
an infinitely large system. However, for systems of finite 
size this threshold p mc (L) depends on the system size. 
The correlation function g(r) for the percolation prob- 
lem is defined as the probability that a site at a distance 
r apart from an occupied site belongs to the same cluster. 
For p < p mc the correlation function is expected to decay 
exponentially as g{r) ~ exp(— r/£) where the correlation 
length £, a measure of the typical cluster diameter, di- 
verges as £ ~ (p mc — p)~ v where, v is the correlation 
length exponent for the diagenetic percolation. 

We use the standard method of estimating the value 
of the percolation threshold. Using a specific sequence of 
random numbers, the lattice is filled at some high value 
of p = phi such that its SC has an infinite cluster. Sim- 
ilarly using the same sequence of random numbers, the 
lattice is filled at some low value of p = p\ so that its 
corresponding SC does not have an infinite cluster. It 
is then similarly tried at a p — (pm + pi )/2. If its SC 
is connecting then phi is equated to p, otherwise pi is 
equated to p. This process is continued till the difference 
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(Phi — Pio) is less than a certain pre-assigned small num- 
ber e = 10~ 5 when p(seq) = (phi +pi )/2 is taken for 
the percolation threshold for this particular sequence of 
random numbers. Averaging over the p(seq) values for a 
large number of independent random number sequences 
one obtains the estimate for p mc (L). 
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FIG. 7. The plot of the percolation probability P2(p,L) 
for m = 2 of the diagenetic percolation for three different sys- 
tem sizes L = 256 (circle), 512 (square) and 1024 (triangle). 



In this process, we tune the probability p to the per- 
colation threshold p mc (L) on a system of size L so that 
the correlation length is of the same order as the system 
size. Therefore for m = 2, L ~ (p 2c — P2c{L))~ u which 
implies 

P2c(L) =p 2c + A.L- 1 ' 1 ' (4) 

We plot p2c{L) in Fig. 4 with L~ l / V and we try v = 4/3, 
the value for the correlation length exponent in the ordi- 
nary percolation. We observe a linear variation for large 
L values. On extrapolation, we find a slightly larger value 
of P2c — 0.5005(2) for sequential updating of type (I) as 
stated above. The p m c{L) values for L=2048, 2896 and 
4096 are found larger than 1/2. For the random and sub- 
lattice sequential updatings the p2 C values are 0.5013 and 
0.5009 respectively. These values should compared to the 
ordinary percolation threshold of 0.592746 on square lat- 
tice 

Since m = 2 is the middle point of the five possible 
values of the neighbour numbers on a square lattice (i.e. 
from to 4) and due to the equivalence of vacant and 
occupied sites, it may be expected that p 2c should be ex- 
actly equal to 1/2. However, we argue that the value of 
P2c very close to 1/2 is actually accidental and there is 
no reason why it should be 1/2. We believe that first ap- 
pearnce of the global connectivity through occupied sites 



determining the percolation threshold is a very special 
situation and since we want this connectivity through 
the occupied sites we break the symmetry between the 
occupied and vacant sites. An estimate of the similar 
percolation threshold for m = 3 on a simple cubic lat- 
tice gives a value much lower than 1/2. However, for 
m = 3 on the triangular lattice the value of the percola- 
tion threshold is obtained as 0.50 ± 0.01. 




[p-p 2c (L)]L- 

FIG. 8. The diagenetic percolation probability P2(p,L) 
as shown in the previous figure for different system sizes are 
scaled as P2(p, V)LP' V with [p — p2 C (L)]L 1 . The collapse is 
obtained using i/=4/3 and (3/v =0.125. 



The p(seq) values corresponding to different sequences 
of random numbers are spread around their mean value 
Pmc{L). The root mean square deviation from the aver- 
age value 

A(L) = (< p(seq) 2 > -[p mc (L)} 2 ) 1/2 (5) 

is supposed to have a dependence on the system size L 
as: A(L) ~ L~ x l v . Plotting A(L) vs. L for m = 2 on a 
double logarithmic scale which fits nicely to a straight 
line gives a value for the correlation length exponent 
v = 1.35(2). 

The fractal dimension df of the "infinite" incipient 
cluster (IIC) of the SC exactly at the percolation thresh- 
old is also calculated. A large number of SCs are gen- 
erated at p = p2c- The average size of the IIC is 
calculated in two ways: (i) Average size S^L) of the 
infinite clusters is measured over the spanning SCs only 
(ii) Average size (L) of the largest cluster is calculated 
over all SCs. Using the definition of percolation proba- 
bility as defined below, S%,(L) = L 2 P 2 (p2 C {L), L). Both 
measures of the IIC are expected to give the fractal di- 
mension: SU 2 (L) - L d f. In Fig. 5 we plot both S^L) 
and 2S 2 C (L) with L on a double logarithmic scale for the 
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system sizes varying from 16 to 4096. The average slopes 
are 1.863 and 1.860 for S^L) and 5^(L) respectively. 
Further, we plot the local slopes d/(L) with 1/L in the 
inset of Fig. 5. After considerable variation over the 
small systems the fractal dimension seems to converge at 
1.89 ± 0.02 for the large system sizes compared to 91/48 
of the ordinary percolation pj] . 




FIG. 9. A stable configuration (SC) at the percolation 
threshold for a system size of L = 80 and with m=3. Sites on 
the "infinite" incipient cluster are joined by lines. 



The cluster size distribution of occupied sites on the 
SCs are also measured at the percolation threshold. We 
define Prob2(S', L) as the probability of a cluster of 
S occupied sites on a SC of a system of size L with 
772 = 2. We start from many independent configurations 
at p2c{L) f» 0.5001 for L = 2048. These are sub-critical 
configurations for the ordinary percolation. We measure 
the Prob2(S l , L) at each time step and keep track of how 
this distribution changes from the initial exponential dis- 
tribution to the power law distribution as shown in Fig. 
6. We notice that at very short times of the order of 1, the 
distribution takes the form of the steady state distribu- 
tion. In this distribution we do not include the "infinite" 
cluster spanning the system. As expected the distribu- 
tion appears to be a power law: Prob2(S l , L) ~ S~ T where 
t = 2.02 ± 0.06 is obtained compared to 187/91 for the 
ordinary percolation pi . 



1.00 r- 

0.98 

0.96 
^ 0.94 
£ 0.92 

0.90 

0.88 




0.86 1 1 1 1 1 1 1 1 1 1 1 

0.00 0.05 0.10 0.15 0.20 0.25 0.30 

Mog(L) 

FIG. 10. Plot of p3 C (L) with l/log(L) which on extrapo- 
lation to L — > oo gives P3 C = 0.96 ± 0.01. 



The order parameter is the percolation probability 
Pm (p) that is the average fraction of sites on the largest 
occupied cluster in the SC. For finite systems it is denoted 
by P m (jp, L) . Variation of the percolation probability is 
shown in Fig. 7 and it varies as: 



P m (p,L) ~ [p — p mc (L)Y 



(6) 



This variation is true in the limit of L — > oo. For fi- 
nite systems however, according to the scaling theory [Q, 
the scaling variable should be where the correla- 

tion length is defined as £ = [p — p m c{L)]~ w . Therefore 
for a finite system of size L the variation of percolation 
probability should be: 



P m (p, L) = L~^ v F[(p - ^(L))! 1 /"] 



(7) 



where, the scaling function F(x) — > ar for large L. We 
show the collapse of the data in Fig. 8 using this scal- 
ing formulation. We again try v = 4/3 and then ob- 
tain a value of f3/v — 0.125 for the data collapse, giving 
j3 = 0.166 compared to 5/36 for the ordinary percolation 

0- 

Next we studied the case of m=3 on the square lat- 
tice. In this case the SC can only be completely vacant 
or it can have only one infinite cluster but cannot have 
isolated clusters. Since in general there will always be 
some sites which have less than 3 occupied neighbours 
on the surface of an isolated cluster, these sites will be 
unstable under the diagenetic rules and the cluster will 
therefore cannot survive in SC. In Fig. 9 we show the 
picture of a SC for 777=3. It is a simple spanning cluster 
having many rectangular holes as in BPM [jl6| . The per- 
colation threshold P3 C (L) also has L dependence and on 
extrapolation with l/log(Z) (as was done in BPM) we 
get p 3 c = 0.96 ± 0.01 (Fig. 10). 
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IV. SUMMARY 

The restructuring process of diagenesis in sedimentary 
rocks involving cementation and dissolution has been 
studied by a percolation model. Simulations on a square 
lattice shows that the porosity is highly reduced due to 
restructuring as is observed in rocks. We also observe 
that starting from the sub-critical configurations of ordi- 
nary percolation at a certain threshold value pi c of the 
pore probability the system evolves to a globally con- 
nected porous space at the stable state. This configura- 
tion is critical since it shows long range correlations. Our 
numerical results give strong indications that the stable 
states in this model have the same critical behaviour as 
that of ordinary percolation. We view the dynamics un- 
der diagenetic rules as a self-organizing dynamics in a 
limited sense since one has to tune p to arrive at a spe- 
cific sub-critical configuration at p2 C so that it organizes 
to show criticality in the stable state. 

After completion of a draft of this manuscript we were 
informed by the editor that a spin model on a square 
lattice where each spin ±1 was flipped only when more 
than half of its four neighbours point into the opposite 
direction was studied in (2l| . Using a much bigger system 
size (Lk7x 10 5 ) compared to what we used a percola- 
tion threshold of 0.5007 ± 0.0001 was estimated which is 
consistent with our results. 
Electronic Address: manna@boson.bose. res. in 
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